#分类别，画252个位点的beta0展示图。
library(ggplot2)
library(dplyr)
library(multiplot)

setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/bayes_pvalue_beta0/")
rt=read.table("../20201110找到有差异的且具有相同pattern的位点/DC.BF1.intersect.txt",head=T,sep="\t")
rt$num=rowSums(rt[,seq(5,25,4)]>=0,na.rm = TRUE)
rt=rt[rt$num>1,]
id=rt$unitID

group1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49","M28_M27","M30_M29","M26_M25","M35_M36","M18_M17","M20_M19","M22_M21","M40_M39")
i=1
fn1=paste0(group1[i],".bayes_p.txt")
file=read.table(fn1,head=T,sep = "\t")
file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
file=file[file$normal_bayes_pvalue<0.05|file$tumor_bayes_pvalue<0.05,]
file=file[file$unitID %in% id,]
rt=file
names(rt)=c("unitID",paste(group1[i],names(file)[2:5],sep="_"))

for (i in 2:6) {
  fn1=paste0(group1[i],".bayes_p.txt")
  file=read.table(fn1,head=T,sep = "\t")
  file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
  file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
  file=file[file$normal_bayes_pvalue<0.05|file$tumor_bayes_pvalue<0.05,]
  ftmp=file[file$unitID %in% id,]
  names(ftmp)=c("unitID",paste(group1[i],names(file)[2:5],sep="_"))
  rt=merge(rt,ftmp,by="unitID",all=T)
}
file=rt
j=ncol(file)-1
file$mean.normal.beta0=rowMeans(file[,seq(2,j,4)],na.rm = T)
file$mean.tumor.beta0=rowMeans(file[,seq(4,j,4)],na.rm = T)
std <- function(x) {
if(length(x[!is.na(x)])==1){return (0)
}else(return(sd(x,na.rm = T)/sqrt(length(x[!is.na(x)]))))
}###算标准误差
file$std.normal.beta0=apply(file[,seq(2,j,4)], 1, std)
file$std.tumor.beta0=apply(file[,seq(4,j,4)],1,std)
file$normal.beta0.upper=file$mean.normal.beta0+file$std.normal.beta0
file$nomal.beta0.lower=file$mean.normal.beta0-file$std.normal.beta0
file$tumor.beta0.upper=file$mean.tumor.beta0+file$std.tumor.beta0
file$tumor.beta0.lower=file$mean.tumor.beta0-file$std.tumor.beta0

#test=tidyr::separate(rt,unitID,into=c("chr","pos","ref","alt"),sep=":")
#test=data.frame(test[,1:2],test[,2:4])
#write.table(test,"../20201110找到有差异的且具有相同pattern的位点/repeat.more.than.2.txt",quote=F,row.names=F,sep="\t")
test=read.table("../20201110找到有差异的且具有相同pattern的位点/DC.BF1.intersect.txt",head=T,sep="\t")
test$num=rowSums(test[,seq(5,25,4)]>=0,na.rm = TRUE)
test=test[test$num>1,]
test=data.frame(unitID=test$unitID,group=test$pattern,repeatnum=test$num)
file=merge(file,test,by="unitID")
anno=read.csv("../20201110找到有差异的且具有相同pattern的位点/repeat.more.than.2.anno.hg19_multianno.csv",head=T)
anno$unitID=paste(anno$Chr,anno$Start,anno$Ref,anno$Alt,sep=":")
anno=data.frame(unitID=anno$unitID,anno[,6:8])
file=merge(anno,file,by="unitID")
rt1=data.frame(file$unitID,file$avsnp150,file$Gene.refGene,file$group,file$mean.normal.beta0,file$normal.beta0.upper,file$nomal.beta0.lower,file$mean.tumor.beta0,file$tumor.beta0.upper,file$tumor.beta0.lower)
names(rt1)=c("unitID","avsnp150","Gene.refGene","group",paste0("DC.",c("mean.normal.beta0","normal.beta0.upper","nomal.beta0.lower","mean.tumor.beta0","tumor.beta0.upper","tumor.beta0.lower")))

normaldata=select(file,unitID,mean.normal.beta0,normal.beta0.upper,nomal.beta0.lower,avsnp150,group)
normaldata=arrange(normaldata,group,mean.normal.beta0)
normaldata$num=1:dim(normaldata)[1]
orders=select(normaldata,unitID,group,num)
normaldata$na=""
p1=ggplot(normaldata,aes(y=mean.normal.beta0,x=num))+scale_x_continuous(breaks = normaldata$num,labels = normaldata$na,position="top")+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=nomal.beta0.lower,ymax=normal.beta0.upper,color=group),width=0.5)+labs(x="",y="con β0")+
  geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c("#1B9E77","#D95F02","#7570B3","#E7298A","#66A61E","#E6AB02","#A6761D","#666666"))+coord_flip()+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line =element_line(size = 0.5))+theme(legend.position="none")
p1
tumordata=select(file,unitID,mean.tumor.beta0,tumor.beta0.upper,tumor.beta0.lower,avsnp150,group)
tumordata$na=""
tumordata=merge(tumordata,orders,by="unitID")

p2=ggplot(tumordata,aes(y=mean.tumor.beta0,x=num))+scale_x_continuous(breaks = tumordata$num,labels = tumordata$na)+
    geom_point(aes(color=group.x),size=1.5)+geom_errorbar(aes(ymin=tumor.beta0.lower,ymax=tumor.beta0.upper,color=group.x),width=0.5)+
    labs(x="",y="case β0")+geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c("#1B9E77","#D95F02","#7570B3","#E7298A","#66A61E","#E6AB02","#A6761D","#666666"))+theme_classic(base_size = 15)+
    theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line = element_line(size = 0.5))+coord_flip()+theme(legend.position="none")
p2

#write.csv(normaldata,"../20201110找到有差异的且具有相同pattern的位点/DC.normal.csv",quote=F,row.names=F)
#write.csv(tumordata,"../20201110找到有差异的且具有相同pattern的位点/DC.tumor.csv",quote=F,row.names=F)
###查看252个位点在双发中的情况
i=7
fn1=paste0(group1[i],".bayes_p.txt")
file=read.table(fn1,head=T,sep = "\t")
file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
file=file[file$normal_bayes_pvalue<0.05|file$tumor_bayes_pvalue<0.05,]
file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
file=file[file$unitID %in% id,]
rt=file
names(rt)=c("unitID",paste(group1[i],names(file)[2:5],sep="_"))

for (i in 8:10) {
  fn1=paste0(group1[i],".bayes_p.txt")
  file=read.table(fn1,head=T,sep = "\t")
  file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
  file=file[file$normal_bayes_pvalue<0.05|file$tumor_bayes_pvalue<0.05,]
  file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
  ftmp=file[file$unitID %in% id,]
  names(ftmp)=c("unitID",paste(group1[i],names(file)[2:5],sep="_"))
  rt=merge(rt,ftmp,by="unitID",all=T)
}
file=rt
j=ncol(file)-1
file$mean.normal.beta0=rowMeans(file[,seq(2,j,4)],na.rm = T)
file$mean.tumor.beta0=rowMeans(file[,seq(4,j,4)],na.rm = T)
std <- function(x) {
if(length(x[!is.na(x)])==1){return (0)
}else(return(sd(x,na.rm = T)/sqrt(length(x[!is.na(x)]))))
}###算标准误差
file$std.normal.beta0=apply(file[,seq(2,j,4)], 1, std)
file$std.tumor.beta0=apply(file[,seq(4,j,4)],1,std)
file$normal.beta0.upper=file$mean.normal.beta0+file$std.normal.beta0
file$nomal.beta0.lower=file$mean.normal.beta0-file$std.normal.beta0
file$tumor.beta0.upper=file$mean.tumor.beta0+file$std.tumor.beta0
file$tumor.beta0.lower=file$mean.tumor.beta0-file$std.tumor.beta0

file=merge(file,orders,by="unitID")
file=merge(anno,file,by="unitID")
rt2=data.frame(file$unitID,file$mean.normal.beta0,file$normal.beta0.upper,file$nomal.beta0.lower,file$mean.tumor.beta0,file$tumor.beta0.upper,file$tumor.beta0.lower)
names(rt2)=c("unitID",paste0("CC.",c("mean.normal.beta0","normal.beta0.upper","nomal.beta0.lower","mean.tumor.beta0","tumor.beta0.upper","tumor.beta0.lower")))

normaldata=select(file,unitID,mean.normal.beta0,normal.beta0.upper,nomal.beta0.lower,avsnp150,group,num)
normaldata$na=""
p3=ggplot(normaldata,aes(y=mean.normal.beta0,x=num))+scale_x_continuous(breaks = normaldata$num,labels = normaldata$na,position="top")+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=nomal.beta0.lower,ymax=normal.beta0.upper,color=group),width=0.5)+labs(x="",y="con β0")+
  geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c("#1B9E77","#D95F02","#7570B3","#E7298A","#66A61E","#E6AB02","#A6761D","#666666"))+coord_flip()+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line =element_line(size = 0.5))+theme(legend.position="none")
p3

tumordata=select(file,unitID,mean.tumor.beta0,tumor.beta0.upper,tumor.beta0.lower,avsnp150,group,num)
tumordata$na=""
p4=ggplot(tumordata,aes(y=mean.tumor.beta0,x=num))+scale_x_continuous(breaks = tumordata$num,labels = tumordata$na)+
    geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=tumor.beta0.lower,ymax=tumor.beta0.upper,color=group),width=0.5)+
    labs(x="",y="case β0")+geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c("#1B9E77","#D95F02","#7570B3","#E7298A","#66A61E","#E6AB02","#A6761D","#666666"))+theme_classic(base_size = 15)+
    theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line = element_line(size = 0.5))+coord_flip()+theme(legend.position="none")
p4

#write.csv(normaldata,"../20201110找到有差异的且具有相同pattern的位点/CC.normal.csv",quote=F,row.names=F)
#write.csv(tumordata,"../20201110找到有差异的且具有相同pattern的位点/CC.tumor.csv",quote=F,row.names=F)
###查看健康中的情况
#查看252个位点在双发中的情况
i=11
fn1=paste0(group1[i],".bayes_p.txt")
file=read.table(fn1,head=T,sep = "\t")
file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
file=file[file$normal_bayes_pvalue<0.05|file$tumor_bayes_pvalue<0.05,]
file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
file=file[file$unitID %in% id,]
rt=file
names(rt)=c("unitID",paste(group1[i],names(file)[2:5],sep="_"))

for (i in 12:14) {
  fn1=paste0(group1[i],".bayes_p.txt")
  file=read.table(fn1,head=T,sep = "\t")
  file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
  file=file[file$normal_bayes_pvalue<0.05|file$tumor_bayes_pvalue<0.05,]
  file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
  ftmp=file[file$unitID %in% id,]
  names(ftmp)=c("unitID",paste(group1[i],names(file)[2:5],sep="_"))
  rt=merge(rt,ftmp,by="unitID",all=T)
}
file=rt
j=ncol(file)-1
file$mean.normal.beta0=rowMeans(file[,seq(2,j,4)],na.rm = T)
file$mean.tumor.beta0=rowMeans(file[,seq(4,j,4)],na.rm = T)
std <- function(x) {
if(length(x[!is.na(x)])==1){return (0)
}else(return(sd(x,na.rm = T)/sqrt(length(x[!is.na(x)]))))
}###算标准误差
file$std.normal.beta0=apply(file[,seq(2,j,4)], 1, std)
file$std.tumor.beta0=apply(file[,seq(4,j,4)],1,std)
file$normal.beta0.upper=file$mean.normal.beta0+file$std.normal.beta0
file$nomal.beta0.lower=file$mean.normal.beta0-file$std.normal.beta0
file$tumor.beta0.upper=file$mean.tumor.beta0+file$std.tumor.beta0
file$tumor.beta0.lower=file$mean.tumor.beta0-file$std.tumor.beta0

file=merge(file,orders,by="unitID")
file=merge(anno,file,by="unitID")
rt3=data.frame(file$unitID,file$mean.normal.beta0,file$normal.beta0.upper,file$nomal.beta0.lower,file$mean.tumor.beta0,file$tumor.beta0.upper,file$tumor.beta0.lower)
names(rt3)=c("unitID",paste0("HC.",c("mean.normal.beta0","normal.beta0.upper","nomal.beta0.lower","mean.tumor.beta0","tumor.beta0.upper","tumor.beta0.lower")))

normaldata=select(file,unitID,mean.normal.beta0,normal.beta0.upper,nomal.beta0.lower,avsnp150,group,num)
normaldata$na=""
p5=ggplot(normaldata,aes(y=mean.normal.beta0,x=num))+scale_x_continuous(breaks = normaldata$num,labels = normaldata$na,position="top")+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=nomal.beta0.lower,ymax=normal.beta0.upper,color=group),width=0.5)+labs(x="",y="con β0")+
  geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c("#1B9E77","#D95F02","#7570B3","#E7298A","#66A61E","#E6AB02","#A6761D","#666666"))+coord_flip()+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line =element_line(size = 0.5))+theme(legend.position="none")
p5

tumordata=select(file,unitID,mean.tumor.beta0,tumor.beta0.upper,tumor.beta0.lower,avsnp150,group,num)
tumordata$na=""
p6=ggplot(tumordata,aes(y=mean.tumor.beta0,x=num))+scale_x_continuous(breaks = tumordata$num,labels = tumordata$na)+
    geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=tumor.beta0.lower,ymax=tumor.beta0.upper,color=group),width=0.5)+
    labs(x="",y="case β0")+geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c("#1B9E77","#D95F02","#7570B3","#E7298A","#66A61E","#E6AB02","#A6761D","#666666"))+theme_classic(base_size = 15)+
    theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line = element_line(size = 0.5))+coord_flip()+theme(legend.position="none")
p6

#write.csv(normaldata,"../20201110找到有差异的且具有相同pattern的位点/HC.normal.csv",quote=F,row.names=F)
#write.csv(tumordata,"../20201110找到有差异的且具有相同pattern的位点/HC.tumor.csv",quote=F,row.names=F)
layout <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 1)
multiplot(plotlist=list(p1,p2,p3,p4,p5,p6),layout=layout)